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I.         INTRODUCTION 

Bi-orthogonality  relationships  for  porous  media  have  been  developed  by  Pro- 
fessors Scandrett  and  Frenzen  [Ref.  1].  They  showed  that  Fraser's  [Ref.  2]  bi- 
orthogonality  relationship  for  the  propagation  of  Rayleigh-Lamb  modes  in  a  plate 
with  traction-free  surfaces  is  a  special  case  of  the  porous  case  in  which  the  porous  slab 
has  zero  porosity.  Bi-orthogonality  relationships  can  be  used  to  solve  scattering  prob- 
lems from  interfaces  of  discontinuity  between  vertically  stratified  porous/elastic/fluid 
media.  In  this  thesis,  we  will  examine  the  bi-orthogonality  relationship  for  the  elastic- 
fluid  case. 

Before  the  bi-orthogonality  relationship  can  be  used,  the  eigenvalues  and  eigen- 
functions  corresponding  to  the  boundary  value  problem  must  be  found.  In  this  thesis, 
we  will  attempt  to  find  these  numerically.  The  question  we  ask  ourselves  is  which 
numerical  approach  is  the  best.  One  method  is  to  determine  Rayleigh-Lamb  modes 
[Ref.  3]  in  the  elastic  solid  and  try  to  perturb  those  solutions  to  obtain  modes  for 
the  fluid  loaded  layer.  The  problem  with  this  method  is  that  modes  may  be  missed. 
Another  approach  is  to  use  the  Rayleigh-Ritz  method  [Ref.  4],  which  is  based  upon 
the  calculus  of  variations.  It  uses  an  expansion  of  elementary  functions  to  "minimize" 
a  functional  to  the  corresponding  boundary  value  problem.  The  functions  must  be 
carefully  chosen  due  to  the  complexity  of  the  resulting  coefficient  matrix  which  will 
have  terms  involving  the  unknown  eigenvalues. 

The  method  chosen  in  this  thesis  is  the  finite  difference  method.  This  method 
is  very  useful  when  numerically  modeling  boundary  value  problems  involving  deriva- 
tives. Using  the  finite  difference  approximations,  we  produce  a  discrete  linear  operator 
which  approximates  the  differential  operator  of  the  boundary  value  problem.  We  can 
then  use  MATLAB  [Ref.  5]  to  solve  for  eigenvalues  of  the  linear  operator  which 
are  theoretically  close  to  the  eigenvalues  of  the  differential  operator.  We  will  also 
use  additional  numerical  methods,  to  be  discussed  later,  in  an  attempt  to  refine  our 


approximate  eigenvalues. 

Chapter  II  will  provide  a  background  in  acoustical  theory  and  will  introduce 
terms  that  the  reader  will  need  to  be  familiar  with  for  further  chapters.  Also,  the 
concept  of  bi-orthogonality  is  explained  in  detail  and  a  simple  example  is  presented. 
Finally,  a  description  of  numerical  techniques  and  methods  that  are  used  throughout 
the  thesis  will  be  given. 

In  Chapter  III,  background  on  the  elastic-fluid  problem  is  introduced.  Detailed 
analytical  solutions  to  several  elastic  half  space  problems  are  provided  to  establish  the 
basis  of  the  numerical  solution.  The  fluid  loaded  plate  case  will  be  the  case  studied 
in  the  numerical  portion  of  this  thesis. 

Chapter  IV  examines  the  numerical  difference  approximations  used  in  the 
problem.  The  modal  wave  speeds  are  presented  and  compared  to  the  "refined"  wave 
speeds.  Comparisons  will  also  be  made  between  the  analytical  and  numerical  eigen- 
functions/vectors  of  the  fluid  loaded  elastic  plate. 

Finally,  Chapter  V  will  provide  conclusions  and  discuss  areas  of  further  re- 
search. 


II.         THEORY  AND  METHODS 

A.      ACOUSTIC  THEORY 

1.        Basic  Acoustics  and  the  Wave  Equation 

Acoustics  is  generally  associated  with  sound.  However,  the  definition  of  acous- 
tics is  the  generation,  transmission,  and  reception  of  energy  that  is  in  the  form  of 
vibrational  waves.  Audible  waves  are  between  20  and  20,000  Hz  where  \Hz  = 
Icycle I second.  The  simple  oscillator[Ref.  6],  or  basic  spring  mass  problem,  is  the 
most  common  example  used  to  demonstrate  basic  acoustics.  The  simple  oscillator  is 
governed  by  the  equation 

d2x 

+  u#r  =  0  (II.l) 


dt2 


with  a  general  solution  of 


x  —  Acos(uj0t)  +  Bsin(co0t)  (H-2) 

where  u>o  is  the  angular  frequency  in  radians  per  second  (rads/sec).   Since  there  are 
2n  radians  per  cycle,  we  define  frequency  as 

/-£.  (H.3) 

In  addition,  we  note  that  the  period  T  is  equal  to  — .    If  the  initial  velocity  and 
displacement  are  known,  we  can  obtain  an  exact  solution  to  Equation  (II.l). 

The  fluid-elastic  problem  that  is  the  central  problem  of  this  thesis  is  a  two- 
dimensional  problem,  and  thus  the  fluid  and  elastic  potential  functions  are  governed 
by  the  two-dimensional  wave  equation 


1  <92<S> 


V2*  =  "^  ("-4) 


c 


dt 


We  assume  an  incident  wave  is  travelling  in  the  positive  ar-direction,  which  strikes  an 
interface  generating  reflected  and  transmitted  waves.  A  plane  time  harmonic  wave 
propagating  in  the  postive  rc-direction  has  the  form 

$  =  ip(y)e-t{wt-k>x)  (II.5) 


LO 


where  k\  =  — .  Substituting  Equation(II.5)  into  Equation(II.4),  one  finds  the  ampli- 
tude  function  must  satisfy  the  relation 


9"{y)  +  {k\  -  k2)<p(y)  =  0  (II.6) 

which  is  the  canonical  form  of  the  boundary  value  problem  to  be  solved. 

2.        Stress-Strain  Relationships  For  Elastic  Solids 

Elastic  solids  have  the  property  that  when  external  pressures  are  removed,  the 
solid  returns  to  its  natural  state.  Strain  is  a  measure  of  displacement  while  stress  is 
a  measure  of  force  per  area.  We  define  the  strain  as 


where  u  is  a  displacement  vector  and  tensor  notation  has  been  adopted.  The  term 
U{tj  means  the  partial  derivative  with  respect  to  Xj,  or  for  example,  u1;2  corresponds 
to  — — .  The  strain,  £tJ,  corresponds  to  elements  of  a  matrix.  Note  that  Eij  is  dimen- 
sionless.  Stress,  r2J,  has  the  dimensions  of  force  per  area.  For  small  strain,  e,  the 
stress,  r,  is  proportional  to  strain.  This  relationship  is  known  as  Hooke's  Law  [Ref. 
7],  which  can  be  written  for  an  isotropic  homogeneous  solid  as 


T{j  —  XSkk^ij  +  2/J,£ij  (II.8) 

where  X  and  /i,  known  as  Lame's  elastic  constants,  have  the  dimension  of  pressure, 
and  ${j  is  the  well-known  Kronecker  delta  defined  as  follows: 


8l3  =  { 


1      i  =  i 

(11.9) 

0     i±  j 


The  displacement  vectors  can  be  represented  in  terms  of  potential  functions. 
The  displacement  equation  of  motion  for  an  isotropic  homogeneous  elastic  solid  [Ref. 
7]  can  be  written  as 

^V2u  +  (A  +  /^)VV-u  =  /)u  (11.10) 

We  consider  a  decomposition  of  the  displacement  vector  in  the  form 

u^V^  +  Vx-0  (11.11) 

It  can  be  shown  that  the  displacement  vector  u  satisfies  Equation  (11.10)  if 

1 

■(p 

(11.12) 


vv  =  —  $ 


v2V>  =  — y> 

Crr 

where  c\   =   and  Cy   =   —[Ref.     7].     The  question  of  the  completeness  of 

P  P 

the  solution  in  Equation  (11.11)  is  guaranteed  by  the  following  theorem  given  in 

Achenbach's  book  on  elastic  solids[Ref.  7]. 

Completeness  Theorem:  Let 

ueC2{V  xT) 
f  e  C(V  x  T) 

and  satisfy  the  equation 

//V2u  +  (A  +  /i)VV  •  u  +  pi  =  pii 
in  a  region  of  space  V  and  in  a  closed  time  interval  T.  Also  let 


f  =  c2  VF  +  4V  x  G. 

Then  there  exists  a  scalar  function  ip(x,t)  and  a  vector-valued  function  ip(x,t)  in 
terms  of  which  u(x,  t)  is  represented  by 

u  =  V(^  +  V  x  ^ 
with 

v-v  =  0, 

and  where  ip(x,t)  and  ij)(x,t)  satisfy  the  inhomogeneous  wave  equations 

V2<p  +  F  =  —(p 

CL 

v2v>  +  G  =  4rV> 

The  proof  of  this  theorem  can  be  found  in  [Ref.  7]. 

Representation  of  the  displacement  vector  u  in  terms  of  potentials  will  become 
important  in  the  remainder  of  this  thesis.  We  restrict  our  analysis  by  assuming  the 
strain  in  the  z  direction  is  zero,  which  corresponds  to  the  two  dimensional  plane  strain 
case.  This  allows  us  to  specialize  our  decomposition  into  the  following 

y  v  (11.13) 

where  the  vector  valued  displacement  potential  is  reduced  to  a  single  nonzero  scalar 
component.  Our  stress  equation  follows  from  Equation (II. 8)  as 


Txy     =     fJ,(Uy     +     Vx) 

Tyy     =      X(UX     +    Vy)     +    2  fjiVy 

where,  as  before,  the  potentials  ip  and  xp  satisfy  Equations  (11.12) 


(11.14) 


For  an  ideal  fluid  the  stress  in  any  direction  equals  the  pressure  stress  of  the 
fluid,  or  in  equation  form 

TXX     =     Tyy     =     ~P  (11.15) 

B.      BI-ORTHOGONALITY 

We  say  that  two  functions  ^p(x)  and  ip(x)  are  orthogonal  over  an  interval 

0  <  x  <  L,  if   /    ip(x)il){x)dx  =  0.    This  is  analogous  to  the  dot  product  of  two 
Jo 

orthogonal  vectors  equaling  zero  .   The  orthogonality  relationships [Ref.  8]  for  sines 
and  cosines  are  as  follows 


fL    .    mrx    .    mirx  0      n^m 

I     sin sin ax  =  <  (11.16) 

Jo  L  L  l     „  _  m 


2 


fL        nisx       m-KX 
I     cos — —cos — - — ax  =  < 
Jo  L  L 


n  7^  m 
n  =  rn^0  (11.17) 


L 

2 

L     n  =  m  =  0 


rL 


[^    .    mrx       mirx 
/    sin — —cos — - — ax  =  0  (11.18) 

Jo  L  L 

This  orthogonality  relationship  allows  us  to  determine  values  of  coefficients  when 
arbitrary  functions  are  expanded  in  terms  of  these  functions.  For  example,  if  f(x) 
can  be  written  in  the  form 


f(x)  =  £  Bnsin—  (11.19) 

n=l  L 

then  we  can  apply  the  orthogonality  of  sines  over  the  interval  [0,  L]  to  obtain 

fL  mirx  S2^  rL        mrx    .    rmrx 

/     j(x)sin — - — ax  —  >    Bn  \    sin — —  sin — - — ax  (11.20) 

^o  L  ~.        Jo  L  L 

which  leads  to 


/    f{x)si 
_  Jo 


nnx 
in — —dx 


Bn  =  J^-n L —  (II.21) 


J 

Jo 


sin 


WKX 


dx 


or  the  more  familiar  form 


Br 


2    fL 


n  =  -  (    f(x)sinT^dx  (11.22) 

L  Jo  L 

The  analysis  is  similar  for  cosines.  Note  that  we  have  formally  interchanged  summa- 
tion and  integration  in  the  above  steps  which  is  valid  provided  the  series  and  integrals 
converge. 

Bi-orthogonality  relationships  generalize  the  orthogonal  relationships  of  clas- 
sical eigenfunction  expansions.  In  [Ref.  1],  these  relationships  were  derived  for  waves 
propagating  at  a  free  surface  and  along  interfaces  between  porous  and  elastic  media. 
The  general  bi-orthogonality  relationship  for  a  two  dimensional  porous  layer  which  is 
fluid  loaded  is: 

f_h  [t*xub  -  r%vA  +  sAUB]vorousdy  +  £*  [-pAUB]flmddy  =  0  (11.23) 

where  the  fluid  layer  and  porous  layer  thicknesses  are  hi  and  h2  respectively,  and 
s  —  —{3p  where  f3  denotes  the  porosity  and  p  denotes  the  interstitial  fluid  pressure. 
Note  also  that  A  and  B  denote  different  modes  within  the  same  solid.  When  the 
porosity  of  the  elastic  media  is  zero,  the  equation  reduces  to  the  fluid  loaded  elastic 
solid  bi-orthogonality  relationship. 

To  illustrate  how  the  bi-orthogonality  property  can  be  applied,  we  will  examine 
the  simple  case  of  two  fluids,  each  with  different  densities  in  contact  at  an  interface 
along  the  y  axis.  In  this  case,  the  bi-orthogonality  relationship  reduces  to  standard 
orthogonality  between  sines  and  cosines.  A  diagram  is  shown  in  Figure(l), 


V 

. 

I 

Py=0 

IT 
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Fluid  2 

►- 
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-+ 
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Figure  1.  Bi-orthogonality  of  Two  Fluid  Layers 

where  fluid  1  satisfies  homogeneous  Neumann[Ref.  8]  conditions,  py  =  0,  and  fluid  2 
satisfies  homogeneous  Dirichlet  conditions,  p  =  0,  along  the  lateral  sides  of  the  layers. 
From  [Ref.  8]  we  know  the  eigenfunctions  are  sines  for  the  Dirichlet  conditions  and 
cosines  for  the  Neumann  conditions.  We  equate  the  pressure  and  flow  normal  to  the 
interface  to  obtain 


P 
U 


=  £ 


=  £ 


(In  -  Rn)^U^(y) 
where  p^(y)  and  U^\y)  are  defined  as  follows: 


T 

-L  rr 


V<?)(s/) 


-i  2 


■  u>p2 


-i  2 


P 

u 


(11.24) 


PnHv)  =  cos(ny)      and    pj^(y)  =  sin(my) 
Un\y)  =  cos(ny)     and     U^(y)  =  sin(my) 
Applying  the  bi-orthogonality  condition  as  follows 


(11.25) 


[u\l)\pl=p2)dy 
J\V[U'=U2}dy 

r pPlU1  =  U*]dy, 

Jo 


(11.26) 


we  obtain  the  following  two  sets  of  equations,  written  in  matrix  form  for  easier  ma- 
nipulation 

j(/ +  r)  =  kf , 

J(T-R)  =  GTf, 

(11.27) 

G(/  +  i?)  =  QT, 

KT(7  -  #)  =  QT. 
We  now  define  Q  =  (Q_1G)  and  K  =  (J_1K),  which,  if  both  formulations  are  to  give 
identical  answers,  results  in  an  identity  QK  =  I.  Solving  for  Q  and  K,  we  have: 


/    cos(py)sin(ny)dy 

Gvn  =  (Q-!G)  =  Jo  

/    $in2{ny)dy 
Jo 


K-nq  =  (J       K)  = 


/     cos(qy)sin(ny)dy 
Jo 

/     cos2(qy)dy 
Jo 


(11.28) 


leading  to 


ypn  — 

IT 


n 


n2  —  (p  —  1); 


K 


2(2  -  5ql)  (  n 


nq 


7T  \n2  —  (q  —  I)2  J 

which  gives  the  matrix  element  {QK,)     as 


where  p  +  n  must  be  even 


where  q  +  n  must  be  even, 


(11.29) 


/   j  ypn^-nq   —    < 


0         if  p  -f  q  is  odd 
8(2  -S„) 


E 


?/' 


^         r["2-(p-l)2][n2-(g-l)2] 


It  can  be  shown  that  the  above  summation  gives 


if  p  -f  q  is  even. 


(11.30) 
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(GIC)pq  =  Spq  (11.31) 

Knowing  these  relationships  and  given  initial  incident  wave  amplitudes,  /,  we  can 
solve  for  the  reflected  and  transmitted  amplitudes  from  either  one  of  the  pairs  of 
vector  equations  given  by  Equations  (11.27). 

C.      NUMERICAL  METHODS 
1.        Newton's  Method 

Newton's  method,  one  of  the  numerical  methods  used  in  this  thesis,  is  a  widely 

used  method  for  solving  nonlinear  equations.  We  will  briefly  review  the  concepts  and 

give  a  short  example  of  this  method.    Newton's  Method  uses  a  line  tangent  to  the 

curve  of  a  function  and  is  based  on  the  linear  approximation  of  the  function.  We  start 

with  an  initial  guess  Xq,  that  we  believe  to  be  close  to  the  root,  and  determine  J{xq). 

Then  we  move  along  the  tangent  line  until  the  point  intersects  the  x-axis.  This  value 

of  x  now  becomes  our  next  iterate.   In  general  terms,  Newton's  Method  is  given  as 

[Ref.  4] 

^n+l  =  xn  — — —  (11.32) 

J  \Xn) 

where  n  =  0, 1,  2,  —  Newton's  Method  converges  rapidly  when  the  initial  guess  is  in 
the  neighborhood  of  the  root.  In  fact,  Newton's  Method  is  quadratically  convergent. 
Here  is  an  example  [Ref.  4]: 

f(x)  =  3x  +  sin  x  —  ex 

f'(x)  =  3  +  cos  x  —  ex 

We  choose  xo  =  0  and  perform  3  iterations  as  follows 

xi  =  x0  —       '  =  .33333 
*2  =  *i  -  t^  =  -36017 

/  yi) 

xs  =  x2-4tP\  =  .360421 7 
/O2) 
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which  is  accurate  to  7  significant  digits.  Although  the  convergence  is  rapid,  there  are 
other  considerations  that  need  to  be  addressed  with  this  method.  Newton's  Method 
requires  2  function  evaluations  per  step.  However,  the  total  number  of  function 
evaluations  is  generally  the  same  or  less  than  other  methods  and  Newton's  Method 
usually  takes  fewer  iterations  to  determine  a  root  with  high  accuracy.  The  main 
drawback  to  Newton's  Method  is  that  it  will  not  converge  for  certain  functions.  This 
happens  when  the  tangent  line  is  approximately  horizontal  or  when  we  repeat  a  value 
previously  used,  for  example  X\  =  X4.  Although  these  are  strong  considerations, 
this  happens  infrequently.  For  our  needs,  we  will  use  this  method  to  refine  our  wave 
speeds  found  in  Chapter  III.  We  will  also  use  this  method  in  refining  our  eigenvalue 
estimates  found  using  the  finite  difference  method  which  is  discussed  in  Chapter  IV. 

2.        Finite  Difference  Method 

The  finite  difference  method  is  used  in  numerically  approximating  derivatives 
and  integrals  of  functions.  We  can  use  this  ability  to  replace  derivatives  with  difference 
quotients  and  thereby  approximate  solutions  to  difficult  differential  equations.  There 
are  three  basic  types  of  difference  quotients;  backward,  central,  and  forward,  with 
central  being  the  most  accurate.  In  certain  instances,  such  as  endpoints,  backward 
or  forward  differences  are  often  used  due  to  the  lack  of  additional  points  beyond 
the  extremity  of  an  interval  necessary  for  accomplishing  the  central  difference.  In 
any  case,  it  is  customary  to  ensure  that  all  difference  operators  have  the  same  or  like 
order  of  accuracy  to  avoid  extra  work  and  increase  the  possibility  of  error  cancellation. 
Difference  operators  can  also  be  applied  more  than  once  in  order  to  obtain  second 
order  and  higher  difference  approximations  to  higher  order  derivatives.  For  this  thesis, 
we  will  concern  ourselves  with  only  second  order  differences.  The  central  difference 
operators  [Ref.  4]  are  as  follows: 
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dx 

Ttv=u 

In 

d2x 

Xj_|_|          £"Li     i    X^  —  ^ 

dt2  'I=Ii 

^             ' 

(11.33) 

+  0(h2) 

where  h  —  £\x  is  the  stepsize.  Since  forward  and  backward  differences  are  0(h),  we 
need  to  expand  xl+i  and  xt_i  as  follows  to  obtain  formulas  of  higher  order  for  our 
backward  and  forward  difference  operators: 


xt+1  =  Xi  +  x'{h  +  x'{h2  +  x'l'h?  + 
xt_!  =  x{  -  x'xh  +  x'-h2  -  x'-'h3  + 
Re-arranging  terms  we  obtain 


(11.34) 


(11.35) 


=  a,-+1      Xi  _  h„  _  h,„  +         ,  forward  difference 

x-V  I  £ 

x(  =  : h  —x"{ —x"  +  0(h3)         backward  difference 

h  2  6 

We  will  be  able  to  substitute  for  the  terms  x'J ',  and  a:'"  in  our  final  discretization.  This 

will  be  discussed  in  more  detail  in  Chapter  IV. 


13 


14 


III.         DISPERSION  EQUATIONS  FOR 
ELASTIC  HALF  SPACES 

In  this  chapter  we  will  examine  several  elastic  fluid  problems  with  infinite 
and  finite  layers.  In  each  problem  a  relationship  known  as  the  Dispersion  Equation 
is  found.  A  system  is  dispersive  if  the  phase  velocity  of  a  wave  depends  upon  the 
frequency  of  the  wave.  Solving  this  dispersion  relationship  at  a  prescribed  frequency, 
we  can  determine  the  wavespeeds  for  propagating  modes  and  use  this  for  an  analysis 
of  the  system.  These  modes  are  also  the  necessary  constituents  in  applying  the  bi- 
orthogonality  relationships  being  studied. 

A.      INFINITE  ELASTIC  HALF  SPACE  WITH  A  FREE 
SURFACE 


Free  Space 

„x 


Infinite  Elastic 
Half  Space 


Figure  2.  Infinite  Free-Elastic  Half  Space 

The  case  of  an  infinite  elastic  half  space  with  a  free  surface  is  described  by 
Achenbach  [Ref.  7].  We  consider  the  two  dimensional  plane  strain  case.  We  apply 
stress  free  boundary  conditions  at  the  free  surface: 

TXy      —     fl(Uy     "f    VX)      =     0  (HI.l) 

Tyy      =      X(UX     +     Vy)     +    2flVy      =     0  (III. 2) 
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Dividing  Equation  (III. 2)  by  p  and  substituting  our  longitudinal  and  transverse  wave 
speeds,  c2L  and  Cj  we  obtain 

c2Lvy+{c2L-2c2T)ux  =  0  (III.3) 

Assuming  a  solution  for  our  vector  potentials,  ip  and  xp,  as  follows, 


ip  =  <p(y)e-t{wt-kx) 
xp  =  $(y)e-ilwt-kx) 
our  solutions  for  u  and  v  become 

U    =    <P>X   +  Ipy    —    ik(p  +   Ipy 

v  =  <py  -  xpx  =  ipy  -  ikip 
and  the  two  boundary  conditions  can  be  written: 


(III.4) 


(III.5) 


uy  -f  ikv  =  0 
{c2L  —  2Cr)iku  +  c2Lvy  =  0 
Next,  our  potential  functions  must  solve  Helmholtz  equations  in  the  elastic  solid, 


(III.6) 


VV  +  k2L<p  =  0 

Y       LY  (III.7) 

V2ip  +  k$ip  =  0 

Dropping  the  tilde  notation  on  the  "amplitude"functions  of  y  and  the  e~l^x~  x>  de- 
pendence, these  functions  satisfy  the  wave  equation  if 

<p»  -  tfip  =  o, 

'  Y  (III.8) 

ip"  -  ftp  =  0, 

where  rj2  =  k2  —  k\  and  £2  =  k2  —  k\.  This  leads  to  exponential  solutions 

<p  =  Ae-™  +  Bevy 

(111.9) 
xp  =  Ce'ty  +  Dt*y 
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where  we  stipulate  that  Re(£,rj)  >  0.    To  ensure  boundedness  of  the  solutions  as 
y  — »  —  oo,  we  neglect  solutions  that  "blow  up"  to  obtain 


y  =  Bevy 
V>  =  Deiy 
Substituting  into  the  boundary  conditions  we  have 

ikrjLp  +  £2xp  +  ik(r](p  —  ikift)  =  0 
(c\  -  2<%)ik(ikip  +  £V0  +  <&(y2(P  -  ifcZ^)  =  ° 
Evaluating  the  amplitude  functions  at  y  =  0  and  simplifying  leads  to 


(111.10) 


(III.ll) 


(111.12) 


{2k2  -  kl)B  -  2ik£D  =  0 

2ikr)B  +  (2k2  -  k%)D  =  0 
The  determinant  of  the  coefficient  matrix  in  Equation  (III.  12)  leads  to  the  Dispersion 
Equation 

(2k2  -  k\)2  -  4k2^V  =  0  (111.13) 

In  this  instance,  the  waves  are  nondispersive  because  the  phase  speed  determined 

by  the  above  relationship  is  independent  of  the  frequency.    The  frequency  can  be 

scaled  out  by  noting  k  =  — ,  and  dividing  the  equation  by  ui  .   Using  ci  =  5690m/3 

c 

and  cj  =  3145m/s,  the  compressional  and  shear  wave  speeds  for  steel,  we  can  apply 
Newton's  Method  to  solve  for  c#,  the  Rayleigh  wave  speed,  and  find 

cR  =  2907m/s 

where  cr  <  cj  <  cl-   From  [Ref.   7]  we  have  the  following  approximate  formula  for 
computing  the  Rayleigh  wave  speed. 

■862  +  1.141/  mT1.. 

cr  w — — cT  (111.14) 
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Using  this  formula,  we  predict  a  Rayleigh  wave  speed  of  2902.25m/s  which  is  very 
close  to  our  result.  Now  letting  D  =  1  we  find  that  B  =  .6748z.  At  the  frequency 
u)  =  1  and  setting  x  =  y  =  0,  the  real  parts  of  the  particle  displacements  (u,  v) 
have  been  plotted  in  Figure(3).  Note  that  this  is  the  time  history  of  the  particle 
displacement  at  the  surface  of  the  elastic  half  space.  In  addition,  as  depth  increases, 
the  amplitude  of  the  Rayleigh  wave  decreases  and  the  motion  alternates  between 
retrograde  and  prograde.  Rayleigh  waves  for  various  elastic  solids  appear  in  Table  (I) 
later  in  this  chapter. 


Rayleigh  Surface 
Waves 


Figure  3.  Rayleigh  Wave  Propagation 


B.      FLUID  LOADED  INFINITE  ELASTIC  HALF  SPACE 


*v 


Infinite  Fluid 
Half  Space 


-x 


Infinite  Elastic 
Half  Space 


Figure  4.  Infinite  Fluid-Elastic  Half  Space 

When  a  fluid  is  added,  we  have  slightly  different  boundary  conditions.  Equa- 
tion (III.  1 )  still  holds,  but  Equation  (III. 2)  is  modified  due  to  the  fluid  pressure  at 
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the  boundary.  This  leads  to 

Tyy      =     \(UX     +     Vy)     +     2  fiVy      =      ~P  (III.  15) 

for  which 

P=-Pf<f>t  (IH.16) 

and  where  (f>  is  the  acoustic  velocity  potential  for  the  fluid  and  pj  is  the  fluid  density. 
The  normal  stress  boundary  condition  becomes 

clvy  +  [c\  -  2c2T)ux  +  iu^  4>  =  0  (111.17) 

Ps 

We  allow  slippage  along  the  fluid/solid  interface  (y  =  0),  but  demand  continuity  of 
the  vertical  velocity  there,  i.e.  vsoud  —  Vfiuld.  This  leads  to  a  third  boundary  condition 

—  iuxfy  —  ukip  +  (f)y  —  0  (III. 18) 

The  velocity  potential  function  must  also  satisfy  a  Helmholtz  equation 

\72<j>  +  k2F<j>  =  0  (111.19) 


u> 


where  kp  =  —  and  cp  is  the  fluid  wave  speed.   As  before  we  determine  the  form  of 
cF 

our  solution  using  the  wave  equation  and  assuming  a  potential  of  the  form 


$  =  ^(y)e-^-^)  (111.20) 

Dropping  again  the  tilde  notation  and  the  e~^wt~  x'  dependence,  the  amplitude  func- 
tion must  satisfy 

<f>yy-  (2(f)  =  0  (111.21) 

where  (2  =  k2  —  kF.  The  amplitude  function  has  the  exponential  solution 
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4>  =  Ee-<y  +  Fe<y  (111.22) 

Disregarding  solutions  that  "blow  up"  as  y  — >  ±00  we  obtain  the  amplitude  functions 
for  the  solid  and  fluid 

<p  =  Bevy 

rp  =  De^  (111.23) 

<f>  =  Ee~Cy 
Substituting  these  into  the  three  boundary  conditions  we  have  at  y  =  0: 

2ikr}(p  +  (£2  +  k2)xp  =  0 
(4 ry2  -  k2{c2L  -  2c2T))ip  -  2ik{c2Trp  +  iio^<f>  =  0  (111.24) 

iuTjif  +  ujkxp  —  (<p  —  0 
Simplifying  we  obtain 

2ikrjB  +  (2k2  -  k2T)D  =  0 
(2fc24  -  k2Lc2L)B  -  2ik£4D  +  tu^E  =  0  (111.25) 

iu7]B  +  ukD  -(E  =  0 
Computing  the  determinant  of  the  matrix  and  simplifying  results  in  the  Dispersion 
Equation.  However,  unlike  the  previous  case,  we  now  have  an  added  term  due  to  the 
presence  of  the  fluid  half  space.  This  equation  is 

{{2k2  -  k2T)2  -  Ak2^]  +  14^  =  0  (111.26) 

C  Ps 

Note  that  once  again  the  frequency  can  be  factored  out  of  the  above  expression  which 
implies  a  nondispersive  media.  Applying  Newton's  Method  using  the  wave  speeds  for 
steel  we  determine  what  is  referred  to  as  the  Scholte[Ref.  7]  wave  speed, 

cs  =  2903.567  -  35.325z     m/s 
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Note  that  in  the  case  with  the  fluid  half  space  we  have  a  complex  wave  speed.  Letting 
£  =  lwe  find  B  =  2.51  +  .395*  and  D  =  -.030  +  3.89z.  Setting  u  =  1  and  plotting 
the  real  part  of  the  horizontal  and  vertical  displacements  at  x  =  y  —  0,  we  obtain  the 
particle  displacement  history  as  this  Scholte  wave  passes  through  the  point  x  =  y  =  0. 
The  plot  looks  the  same  as  Figure  (3). 

C.      FINITE  FLUID  INFINITE  ELASTIC  HALF  SPACE 


Fluid 


^x 


Infinite  Elastic 
Half  Space 


Figure  5.  Finite  Fluid  Infinite  Elastic  Half  Space 

We  next  consider  a  fluid  layer  of  finite  thickness  overlying  an  elastic  half  space. 
The  boundary  conditions  from  before,  Equations  (III. 15),  (III. 17),  and  (III. 18)  still 
hold.  However,  an  additional  boundary  condition  is  needed  at  the  boundary  between 
the  fluid  and  free  space  which  is  that  the  pressure  be  zero  there.  From  Equation 
(III. 16),  this  leads  to 

p  =  iupf(f>  =  0  =>  <f>  =  0  (111.27) 

Also,  from  before  we  obtain  the  solutions  to  the  Helmholtz  equation  where 

¥  =  Aevy 

(111.28) 
0  =  Beiy 

Our  third  potential  changes  because  we  now  have  a  finite  layer.    We  are  unable  to 

discard  one  of  the  solutions,  but  instead  must  have 
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<f>  =  Ce<y  +  De'^  (111.29) 

Simplifying  this  by  combining  the  terms,  we  choose 

<j>=  C smh{({h  -  y))  (111.30) 

Notice  that  by  making  this  choice  we  have  incorporated  our  additional  boundary 
condition  so  that  when  y  =  h,  4>  =  0.  Substituting  into  the  boundary  conditions  at 
y  =  0,  we  have 

2ikri<p  +  (£2  +  fc2)i/>  =  0 
(c2Lrj2  -  k2(c\  -  2c2T))p  -  2ifcf4V  +  i"%<t>  =  0  (111.31) 

iurjtp  +  ukip  —  (C  cosh(((h  —  y))  =  0 


or 


2ikrjA  +  {2k2  +  k2T)B  =  0 
{2k2c2T  -u2)A-  2ikic\B  +  iu%  sinh(C^)C  =  0  (111.32) 

iurjA  +  ukB  -  (  cosh((h)C  =  0 
Computing  the  determinant  and  simplifying  leads  us  again  to  a  Dispersion  Equation, 
which  is 


(2k2  -  k2Tf  -  Ak2r,A  +  *4^—  tanh(CA)  =  0  (111.33) 

L  J  C  Ps 

Unlike  the  previous  case,  the  frequency  can  no  longer  be  factored  out.  This  dispersion 
relation  signals  that  we  now  have  a  dispersive  media  which  is  intermediate  between 
the  other  two  limiting  cases.  Notice  that  as  u  — *•  0,  the  "fluid  term"  is  0(u5)  while  the 
solid  terms  are  0(u4),  and  therefore  the  wave  acts  like  a  Rayleigh  wave.  Asw^  oo, 
t  aula  {(h)  ->  1  we  obtain  the  Scholte  limit.  This  makes  sense  physically  because  as 
u  — ►  0,  A  ~  O(-)  becomes  large  and  the  fluid  layer  will  appear  to  become  extremely 
"thin"  relative  to  the  propagating  wavelength.    On  the  other  hand  if  w  — »  oo,  the 
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fluid  will  appear  to  be  "thicker"  in  terms  of  the  number  of  wavelengths  that  can  fit 
onto  the  fluid  layer,  and  we  get  the  other  limit  of  an  infinite  fluid  layer.  As  in  the 
previous  problems  we  solve  for  the  Scholte  interface  wave  speed.  A  height  of  \rn  is 
used  for  the  fluid  layer  in  the  computation.  The  wave  speed  for  steel  is 

c  =  2906-  .001  \i     m/s 

Letting  C  =  1  we  find  that  A  =  .00947  -  .104:rl0-82  and  B  =  .147xl0-8  -  .0057H. 
We  now  set  w  —  1,  h  =  1,  and  x  =  y  =  0.  This  plot  is  the  same  as  our  Rayleigh  wave 
plot  in  Figure  (3),  which  was  expected.  Listed  below  in  Table  (I)  are  Rayleigh  and 
Scholte  wave  speeds  for  various  elastic  solids. 


Solid 

cl 

cT 

cr 

Cq 

Jtn  f  intte  f  lutd 

CSf,n,teflutd 

Steel 

5690 

3145 

2907 

2904  -  35z 

2906  -  .00 It 

Iron  (Cast) 

4175 

2308 

2133 

2127  -  45i 

2133  -  .0009z 

Aluminum 

6242 

3144 

2930 

2904  -  112t 

2930  -  .003e 

Glass  (Pyrex) 

5637 

3297 

3026 

2993  -  153z 

3026  -  .004z 

Rubber  (Hard) 

2117 

864 

814 

694 

814 

Table  I.  Rayleigh  and  Scholte  Wave  Speeds  for  Various  Elastic  Solids 


D.      FLUID  LOADED  ELASTIC  PLATE 
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y  = 


=  o 


Elastic  Soli< 


y*-h± 


Figure  6.  Fluid  Loaded  Elastic  Plate 

We  consider  an  elastic  layer  of  finite  thickness.    We  examine  the  relationhip 
of  a  finite  fluid  media  in  contact  with  a  finite  elastic  medium  or  plate.    As  shown 
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in  Figure  (6),  the  origin  is  assumed  to  go  through  the  middle  of  the  elastic  solid. 
This  will  help  later  to  show  the  symmetric  and  anti-symmetric  modes  as  discussed  by 
Achenbach  [Ref.  7]  in  the  unloaded,  "free  plate"  problem.  Our  boundary  conditions 
are  slightly  different  with  this  case.  We  choose  a  thickness  hi  for  the  elastic  layer  and 
a  thickness  h2  for  the  fluid  layer.  At  y  =  —  —  we  have 


' xy   —  " 

lyy  V 


(111.34) 


Aty  =  +y, 

we  have 

rxy  =  0 
7yy  =  ~P 

Vsolid  =  V  fluid 

(111.35) 


hi 

Finally,  at  y  =  h?  +  — ,  we  want  the  pressure  to  again  be  equal  to  zero.   We  again 

assume  solutions  to  the  wave  equation.  However,  we  now  have  a  finite  fluid  and  elastic 
layer,  and  therefore  must  retain  all  of  the  exponential  terms.  We  write  the  solutions 
in  terms  of  cosh  and  sinh  in  order  to  better  see  the  symmetric  and  anti-symmetric 
portions  of  the  solutions.  We  have, 

<^  =  A  cosh(r]y)  +  B  sinh^t/) 

0  =  C  cosh(£y)  +  D  sinh(£y)  (111.36) 

<t>  =  Esm\i{ah-y)) 

hi 
where  h  =  —  +  /i2,  and  we  have  taken  into  account  the  boundary  condition  that 

p  =  0  at  h. 

Substituting  into  our  boundary  conditions  and  letting  p  =  — ,  we  have 
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2ikrnpy{-%)  +  (2k2  -  k2T)rP{-^-)  =  0 
2tkV^y(^)  +  (2k2  -  k^)xP(^)  =  0 
(2k2  -  JfeJM-£)  -  2»W„(-^-)  =  0  (111-37) 


(2k2  -  k2T)^)  -  2»Wv(f )  +  ^M  =  0 


CT 


iu>wy(!f)  +  wfcV>(^)  -  <f>y(h2)  =  0 
Trying  to  find  and  simplify  the  determinant  is  difficult.  We  try  a  simpler  approach  by 
adding  and  subtracting  like  equations  to  obtain  our  Dispersion  Equation.  Adding  and 
subtracting  the  first  two  equations  and  doing  the  same  with  the  next  two  equations, 
we  have  the  following 

2ikrjB cosh(r)^-)  +  (2k2  -  k^)C cosh(^)  =  0 
2iJb/>lsinh(j7^)  +  (2k2  -  k^)Dsmh(^)  =  0 
(2k2  -  k$)Acosh(r,!f)  -  2ik(D  cosh(^)  +  l-£fEsmh((h2)  =  0  (111.38) 

(2k2  -  k%)Bsinh{r)!f)  -  2ih(C sinh(^)  +  ^Esmh((h2)  =  0 

iuW>y(%)  +  «>fy(%)  -  <j>y(h2)  =  0 

From  the  first  two  equations,  we  solve  for  A  and  B  in  terms  of  C  and  D.  Substituting 
into  the  next  two  equations  we  obtain  a  relationship  of  C  and  D  and  find  E  in  terms 
of  C .  This  is  not  as  difficult  as  it  seems,  but  does  require  accurate  bookkeeping  to 
keep  up  with  the  algebra.  Substituting  all  of  these  into  the  last  equation,  and  after 
much  algebra,  we  obtain  the  Dispersion  Equation  for  a  fluid  loaded  plate. 


[(2A;2  -  k2T)2  -  4k2rjZ  tanh(r/^)  coth(^)][(2P  -  k2T)2  -  4k2^  cothfa^-)  tanh(^)] 
=  -Pfkh^H(h2)[{2k2  _  fc2)2coth(7//li)_4fc2^coth(^i)] 

PsQ 

(111.39) 
Because  the  dispersion  relationships  will  involve  infinitely  many  complex  roots,  due 
to  the  finite  nature  of  the  elastic/fluid  system,  Newton's  method  may  indeed  be  able 
to  find  one  or  several  of  the  roots,  but  will  not  be  efficient  in  finding  the  first  N  roots 
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(ordered  by  magnitude),  simply  because  required  starting  values  are  not  available. 
One  lacks  the  confidence  that  all  of  the  eigenvalues  will  be  found  with  Newton's 
method  without  sufficiently  close  starting  values. 
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IV.         RESEARCH  RESULTS 

A.      DISCRETIZATION  PROCESS 

We  begin  the  numerical  analysis  of  the  fluid  loaded  plate  by  applying  finite 
difference  approximations  to  the  boundary  value  problem.  We  know  that  in  the  fluid 
and  solid,  the  respective  reduced  Helmholtz  equations  must  be  satisfied.  We  define 
(N  +  l)hj  as  the  thickness  of  the  fluid  layer  and  (M  +  l)hs  as  the  thickness  for  the 
elastic  layer,  where  (N  +  1)  is  the  number  of  discrete  points  in  the  fluid  layer  and 
(M  +  1)  as  the  number  of  grid  points  in  the  elastic  layer,  and  hj  and  hs  are  the 
stepsize  for  the  fluid  and  elastic  layers,  respectively.  In  order  to  obtain  our  unknown 
modal  wave  number,  k2,  on  the  right  hand  side  of  the  system  of  equations,  several 
substitutions  are  necessary.  We  define  xj)*  =  x/)x,  and  also  substitute  4>t  =  4>.  Therefore 
in  the  final  analysis  we  will  be  solving  for  normal  modes  of  the  functions  <#;,  ip,  and 
xp*.  Dropping  the  superscripts  and  subscripts  to  avoid  clutter  in  the  discretization, 
we  have 

#  +  k)<f>i  =  kl<j>i        i  =  \..N-\ 

<rf  +  *i<Pi  =  l£<Pi        i  =  l..M-l  (IV.l) 

0f+*J^  =  A£0<        z  =  l..M-l 
The  additional  boundary  conditions  on  the  fluid  loaded  elastic  plate  from  Chapter 
III,  can  be  reduced  after  completing  the  above  substitutions  and  applying  the  finite 
difference  method,  to  the  following: 

4>N  =  0 

2<p"  +  (2k2L  -  k2T)ip  -  2x1)'  =  ^(j) 

-2k2y  +  k2Txl)  +  2</>"  =  0 


*-•*>'  =  & 

2y"  +  (2k\  -  k$)<p  -  2x1)'  =  0 

-2k2<p'  +  k\x\)  +  2x1)"  =  0 
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(IV.2) 


where  i '  =  0  is  the  interface  between  the  fluid  and  the  solid.  In  order  to  incorporate 
these  additional  boundary  conditions  into  Equations  (IV.  1),  we  must  solve  for  what 
we  will  call,  pseudo  nodes.  A  pseudo  node  is  a  point  outside  of  the  fluid  or  elastic 
solid  that  will  be  used  in  the  differencing  of  a  potential.  This  can  be  seen  below  in 
Figure  (7), 


J 

, 

y 

N 

Fluid 

<t> 

-1 
0    . 

M+  1 


Figure  7.  Discretization  of  Fluid  Loaded  Elastic  Plate 

where  <$>\  is  the  pseudo  node  for  the  fluid  and  </?_i,  0_i,<£>m+1i  and  V'm+i  are  the 
pseudo  nodes  for  the  elastic  solid. 

The  equations  for  the  pseudo  nodes  require  some  careful  algebra  and  can  be 
reduced  to  the  following: 

^  =  ^  +  h-hlfUL  fo  +  M^ii!^  _  2hjuHi  +  ±g£[l  +  h2(k2T  -  kl)]<p0  -  hshfk2Lu2<pi  -  ^£<p2 


kthi 


V?_i  =  [1  +  h2s{k2T  -  k2L)]<p0  +  (1  -  fcJ*£)Pi  -  <P2  +  (2/i,  +  ^)^o  -  2/isVi  +  ^to 


hsPf 


0_,  =  (4  +  fcf.7^o  -  3^i  +  [hsk2T  -  ^)^0  +  (£  -  2A£/i,)y>1  -  ^  V2  +  ^<^ 
9?m+i  =  [1  -  h2{k2L  -  fcf.)]^?Af  +  (1  -  h2k2L)<pM-i  -  Pm-2  +  2hsxl>M-i  +  ( 


h*ki 


2hs)il> 


M 


ipM+i  =  (4  -  h2sk\)\l>M  -  ZxpM-i  +  (ir-  hsk^)(fM  +  {1k2Lhs  -  -g-)<pM-i  +  TT^M-2 


(IV.3) 


These  pseudo  nodes  are  substituted  into  our  Helmholtz  equations  from  (IV.  1)  to 
obtain  equations  for  the  boundaries  i  =  0  and  i  —  M .  Our  discretized  boundary 
equations  are  now  entered  into  a  MATLAB  code  given  in  the  Appendix.  Note  that 
for  the  MATLAB  code  we  must  start  our  index  at  1  as  MATLAB  does  not  have  0 
indices.  A  shift  is  required  so  that  the  first  N  places  correspond  to  <J>,  the  next  M  -f  1 
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places  correspond  to  tp,  and  lastly,  the  next  M  +  1  places  correspond  to  x\).  This  allows 
us  to  keep  track  of  each  of  the  three  potentials.  MATLAB  computes  our  eigenvalues 
and  orders  them  in  ascending  magnitude.  The  code  also  calculates  the  normal  and 
shear  stresses  for  use  in  applying  the  bi-orthogonality  condition.  Finally,  various 
graphs  involving  our  discretized  solutions  will  be  plotted  against  analytical  solutions. 
The  analytical  solutions  for  stresses,  pressures,  and  displacements  are  determined 
using  the  coefficients  listed  in  the  last  section  of  Chapter  III. 

B.      NUMERICAL  RESULTS 
1.        Modal  Wave  Numbers 

The  elastic  solid  used  for  the  following  analysis  will  be  steel.  We  begin  with 
an  examination  of  the  eigenvalues  of  our  discretized  matrix.  In  the  discretization, 
uj  =  27r,  N  =  10,  and  M  =  50,  leading  to  a  square  matrix  of  size  112  from  which  112 
eigenvalues  are  determined.  Since  we  are  solving  for  fc2,  we  must  take  the  square  root 
of  the  eigenvalues  and  place  them  in  order.  Determination  of  the  accuracy  of  these 
wave  numbers  will  involve  substitution  into  Equation  (III. 39),  our  fluid  loaded  elastic 
plate  dispersion  relation.  A  plot  of  the  norm  of  the  eigenvalues  is  shown  in  Figure 
(8)  on  the  following  page.  Aside  from  the  last  few  eigenvalues,  we  get  a  nice  ordering 
of  the  eigenvalues. 

We  attempt  to  refine  our  approximate  wave  numbers  using  several  iterations  of 
MAPLE  coded  Newton's  Method  with  the  eigenvalues  as  initial  iterates.  The  results 
are  shown  in  Table  (II).  For  some  of  our  "refined"  wave  speeds,  Newton's  Method 
attempts  to  converge  to  a  root.  Newton's  Method  diverged  for  the  roots  marked  with 
an  asterik.  After  refinement,  several  of  the  roots  show  improvement  in  satisfying  this 
equation,  i.e.  the  magnitude  of  the  residual  of  the  relationship  is  much  closer  to  zero. 
The  fundamental  root  has  a  residual  magnitude  change  from  10~23  to  10~50.  The 
second  and  third  roots  change  almost  as  much.  However,  the  rest  of  the  roots  exhibit 
less  of  a  change  in  magnitude,  indicating  that  the  starting  values  for  the  roots  may 
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Norm  of  Eigenvalues 


0  0.05  0.1  0.15  0.2  0.25  0.3  0.35  0. 

magnitude  of  lambda 

Figure  8.  Eigenvalues  of  Discretized  Matrix 
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not  be  close  enough  to  the  actual  roots  for  Newton's  Method  to  converge. 


n 

k 

Refined  kn 

relative  error 

1 

.001291  -  .0003723* 

.001218 

0.2824 

2 

.001291  +  .0003723i 

.001218 

0.2824 

3 

.008880 

.004932 

.4445 

4 

.006711  -.0080032* 



5 

.006711  +  .0080032* 

6 

.004703  -  .01666e* 



7 

.004703  +  .016662* 



8 

.002550  -  .03182 

.01127-  .02102 

.4355 

9 

.002550  +  .03182 

.01127 +  .02102 

.4355 

10 

.001664  -  .04762 

.0007991  -  .03982 

.1646 

Table  II.  Calculated  versus  Refined  Modal  Wave  Speeds 

Table  (II)  shows  the  first  three  roots  converging  to  a  real  root. The  other  "refined" 
knS  do  not  appear  to  converge  to  any  of  the  k^s  used  as  initial  starting  values.  It 
may  be  necessary  to  search  more  than  the  first  ten  modes  to  help  determine  if  a  shift 
of  the  "refined"  values  to  a  corresponding  kn  is  possible.  For  this  reason,  we  will 
concentrate  on  the  first  few  modes. 

2.        Stress  Potentials 

The  following  graphs  show  the  analytical  versus  the  discretized  solutions  for 
the  fluid  and  solid  potentials  <£,  (p,  and  ip.  Figure  (9),  displays  the  best  results.  The 
analytical  solution  of  (f>  is  almost  identical  with  the  discretized  solution.  Note  that 
the  discretization  code  has  <p\  as  the  first  interior  point  and  not  the  boundary.  This 
accounts  for  the  slight  mismatch  at  the  fluid- free  space  boundary.  Figure  (10)  shows 
a  simliar  pattern  where  our  solution  is  nearly  identical,  but  diverges  slightly  as  we 
near  the  fluid-elastic  boundary.  Figure  (11)  diverges  at  both  ends  but  is  on  a  scale  of 
10~5.  When  viewed  on  a  small  scale,  the  discretized  solution  does  diverge  more  than 
<f>  and  tp.  It  is  obvious  that  there  is  a  problem  at  the  fluid-elastic  boundary  and  from 
Figure  (11),  there  is  a  slight  indication  of  a  problem  at  the  bottom  boundary. 
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Analytical  vs  Discretized  for  Phi  (fluid) 
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Figure  9.  Analytical  versus  Discretized  Solutions  for  <f> 
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Analytical  vs  Discretized  for  Varphi  (solid) 
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Figure  10.  Analytical  versus  Discretized  Solutions  for  ip 
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Analytical  vs  Discretized  for  Psi  (solid) 
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Figure  (12)  illustrates  a  pressure-stress  graph.  Notice  that  at  the  boundary 
of  the  fluid  and  elastic  layer  a  discontinuity  exists.  One  of  our  boundary  conditions 
from  the  fluid  loaded  plate  is 

Tyy  =  ~P  (IV.4) 

This  means  that  our  plot  should  connect  smoothly  at  the  interface.  The  discontinuity 
indicates  a  problem  with  the  discretization  at  the  interface.  This  problem  was  evident 
in  the  previous  plots  and  is  re-enforced  with  this  graph. 

Figure  (13)  plots  the  wave  numbers  of  the  first  three  modes  as  a  function  of 
increasing  frequency.  We  see  that  our  fundamental  mode  is  purely  real  and  our  second 
and  third  modes  are  complex. 

3.        Bi-orthogonality 

The  bi-orthogonality  condition  of  differing  modes  from  Chapter  II  indicates 
that  we  should  expect  a  diagonal  matrix  in  forming  the  bi-orthogonality  product  of  a 
set  of  modes,  where  the  diagonal  element,  xn,  is  the  "norm"  of  each  individual  mode. 

*n  =  f°   [tWuW  -  r^%lasticdy  -  fhl  \p{n)U(%luiddy  (IV.5) 

J-h.2  JO 

The  bi-orthogonality  condition  states  that  the  off-diagonal  elements,  the  product  of 
two  separate  modes,  should  go  to  zero.  Numerically,  we  expect  these  elements  to  tend 
toward  zero.  The  numerical  bi-orthogonality  relationship  for  the  fluid  loaded  plate 
yields  the  following  matrix: 


4.441  x  105  4.561  x  105  5.233  x  103  3.639  x  103  3.836  x  103 

4.561  x  105  4.441  x  105  5.233  x  103  3.836  x  103  3.639  x  103 

5.118  x  105  5.118  x  105  7.046  x  105  2.670  x  106  2.670  x  106            (IV.6) 

2.617  x  106  1.005  x  107  6.139  x  106  1.524  x  107  1.832  x  107 

1.005  x  107  2.617  x  106  6.139  x  106  1.832  x  107  1.524  x  107 
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Pressure  -  Normal  Stress  Mode  1 
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Figure  12.  Pressure  -  Normal  Stress  in  Layer 
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Wave  Propagation  for  Steel 


Third  Mode 


x10 


Imaginary  ~10      0 


0.04 


Real 


Figure  13.  Wave  Mode  Propagation  -  Modes  1-3 
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Notice  in  this  matrix,  most  of  the  off-diagonal  elements  tend  to  exhibit  some  form 
of  reduction.  However,  this  did  not  happen  as  quickly  as  believed.  This  indicates 
that  the  bi-orthogonality  condition  is  very  sensitive  to  small  errors  in  calculating  the 
modes.  A  more  careful  numerical  treatment  might  improve  the  diagonal  form  of  the 
matrix,  but  it  may  be  necessary  to  limit  the  application  of  the  numerical  method 
to  eigenmode  determination  only,  from  which  the  eigenmodes  can  be  refined  and 
ultimately  found  by  analytical  means. 
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V.         CONCLUSIONS  AND  FUTURE 

RESEARCH 

We  can  draw  several  conclusions  on  the  discretization  of  the  fluid  loaded  elastic 
plate.  The  eigenvalues,  or  modal  wave  numbers,  that  resulted  from  the  discretization 
nearly  satisfy  the  Dispersion  Equation  (III. 39).  Using  Newton's  Method,  we  deter- 
mined that  our  initial  values  may  not  be  close  enough  in  some  instances  to  the  actual 
roots  as  is  evidenced  from  the  lack  of  improvement  of  the  residual  calculation  in  the 
"refined''  root.  Therefore,  we  conclude  that  there  are  numerous  other  roots  that  we 
have  not  found. 

The  biggest  disparity  lies  in  the  plot  of  pressure  and  normal  stress  versus  the 
height  of  the  fluid  and  elastic  layer.  As  seen  in  the  discretization  code  in  the  Ap- 
pendix, differencing  of  the  potentials  is  necessary  when  solving  for  the  stresses  and 
displacements.  Therefore,  any  numerical  error  in  the  fluid-elastic  boundary  differenc- 
ing is  compounded.  Recall,  in  theory  the  normal  stress  is  equal  to  minus  the  pressure 
at  this  boundary.  The  discretization  does  not  show  this.  Finally,  we  conclude  that 
our  graph  of  the  first  three  modes  in  Figure  (13),  could  indeed  be  modes  for  this 
problem.  However,  they  are  most  likely  not  the  first  three  propagating  modes. 

The  trouble  with  the  results  is  mostly  due  to  the  following:  a)  there  are  errors 
in  the  discretization  or  b)  the  discretization  is  unable  to  properly  model  the  boundary 
conditions.  Indications  of  the  first  problem  are  evidenced  in  two  ways.  The  first  is  the 
inability  of  the  eigenvalues  to  converge  to  the  correct  analytical  eigenvalues.  Secondly, 
the  plotting  of  the  corresponding  eigenvectors  do  not  satisfy  the  boundary  conditions 
as  required.  The  solution  to  this  problem  is  to  find  the  correct  way  to  difference  the 
boundary  conditions.  Nothing  can  be  done  to  fix  the  second  problem,  except  to  try 
a  different  discretization.  The  code  as  it  stands  does  find  some  roots,  but  until  the 
boundary  conditions  can  be  satisfied,  there  is  no  real  confidence  that  the  roots  are 
comprehensive. 
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Several  new  problems  arise  from  this  research.  The  first  is  to  determine  a 
better  way  of  differencing  the  boundary  of  the  fluid-elastic  layer.  Once  this  problem 
is  corrected,  we  can  solve  for  the  coefficients  of  transmitted  and  reflected  waves  given 
and  incident  wave.  This  is  useful  when  extending  the  elastic  case  to  the  more  complex 
porous  case.  If  the  porous  case  can  be  solved  numerically,  this  will  allow  us  to  create 
ocean  acoustic  models  where  the  transmitted  and  reflected  waves  can  be  determined 
from  known  coefficients.  This  has  potential  applications  in  sonar  and  mapping  of 
the  ocean  floor.  Finally,  a  similiar  problem  involves  modeling  a  fluid  layer  only,  with 
different  sound  velocity  profiles  and  an  acoustically  hard  boundary.  The  discretization 
would  not  involve  the  elastic  boundary  and  may  be  a  better  way  to  approach  the  fluid- 
elastic  case. 
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APPENDIX.  MISCELLANEOUS  DATA 

1.  MATLAB  CODE  FOR  GENERATING  WAVE  SPEEDS 

function    []    =elastic(E,   nu,    rho) 

°/o  ************************************************************** 
*/,   ELASTIC  -  This  function  computes  the  elastic  constants  and 
°/0  wavespeeds  for  various  solids.  Input  values  are  the 

°/0  Young's  Modulus,  Poisson's  ratio,  and  the  density. 

% 

%   Copyright  1997  by  C.  R.  Myers,  III 

°/o  ************************************************************** 

lambda=0; 

mu=0; 

cl=0; 

ct=0; 

E=E*10~10; 

lambda  =  (E*nu)/( (l+nu)*(l-2*nu)) 
mu  =  E/(2*(l+nu)) 

cl  =  round (sqrt(( lambda  +  2*mu)/rho)); 
ct  =  round (sqrt (mu/rho)) ; 

2.  DISCRETIZATION  MATLAB  CODE 

°/o  Fluid-Elasic  Finite  Difference  Matlab  Code 

°/, 

1   Copyright  1998  by  Coley  R.  Myers,  III 

°/o  Initialize 

clear; 

format  long  e; 

w=0; 

E_val=[]; 

%  allows  program  to  run  in  loop  for  obtaining 
c/o  modes  over  a  specified  frequency  range. 
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counter=l ; 

°/,  loop  is  a  counter  and  is  used  to  obtain  graph 
°/0  graph  of  modes  over  a  specified  frequency  range 

for  loop=l : counter 

°/0  Variables  input  for  specific  elastic  solids 

I=sqrt(-1); 

Lf=10; 

Ls=200; 

w=2*pi+(w+.l);  rf=1026;  rs=7700; 

ct=3145;  cl=5690;  cf=1500; 

modes=10; 

mu=rs*ct~2; 

lam=rs*cl~2  -  2*mu; 

hf=l; 

hs=4; 

N=Lf /hf ; 

M=Ls/hs; 

d=N+M+M+2 ; 

A=zeros(d) ; 

kf=w/cf;  kt=w/ct;  kl=w/cl; 

°/0  Ensures  matrices  are  clear  when  loop  >  1 

Vl=[]  ;Phi=[]  ;Vphi=[]  ;Psi=D  ;Psixy=D  ;Vphiy=[]  ; 

V=[];D=[];Pres=[];U=D; 

u=  []  ;  v=  []  ;  Styy=  []  ;  Stxx=  D  ;  Stxy=  D  ; 

pr.styy=[]  ;f  1=  D  ;F1=D  ;Y=D  ;Y1-D  ; 

Vphia=[]  ;Psia=D  ;Phia=[]  ; 

°/.  Fluid  BC  Phi(0)=0  and   1st   equation 

P-lJ 

A(p,l)=-(2/hf^2-kf"2) ; 

A(p,2)=l/hf-2; 
p=p+l; 

for   i=2:N-l 
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A(p,i-l)=l/hf-2; 
A(p,i)=-(2/hJT2-kf-  2); 
A(p,i+l)=l/hf~2; 

p=p+l; 

end; 

'/.  3rd  BC  for  Phi 

A(p,N-l)=2/hf~2; 

A(p,N)=-2/hf~2+kf~2+(hs*w~2*rf)/(mu*hf) ; 
A(p,N+l)=(w~2/(hs*hf))*(l+hs~2*(kt~2-kl~2)); 
A(p,N+2)=-hs*w~2*kl~2/hf ; 
A(p,N+3)=-w"2/(hs*hf); 
A(p,N+M+2)=(kt~2*hs~2*w-2)/(2*hf); 
A(p,N+M+3)=-2*w~2/hf ; 

•/.  1st  BC  for  Varphi 

p=p+l; 

A(p,N+M+2)=(2/hs)+kt~2*hs/2; 

A(p,N+M+3)=-2/hs; 

A(p,N+3)=-l/hs~2; 

A(p,N+2)=(2/hs~2)-kl~2; 

A(p,N+l)=kt~2-l/hs~2; 

A(p,N)=rf /mu; 

P=P+1; 

'/.  2nd  BC  for  Varphi 

for  i=(N+2):(N+M) 

A(p,i-l)=l/hs~2; 
A(p,i)=-(2/hs~2-kr2); 
A(p,i+l)=l/hs~2; 
p=p+l; 

end; 

1   3rd  BC  for  Varphi 
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A(p,N+2*M+2)=(kt~2*hs/2)-2/hs; 

A(p,N+2*M+l)=2/hs; 

A(p,N+M+l)=kt-2-l/hs"2; 

A(p,N+M)=2/hs~2-kl~2; 

A(p,N+M-l)=-l/hs~2; 

p=p+l; 

X  1st  BC  for  Psi 

A(p,N)=rf/(mu*hs); 

A(p,N+l)=kt^2/hs-2/hs"3; 

A(p,N+2)=4/hs~3-2*kl"2/hs; 

A(p,N+3)=-2/hs~3; 

A(p,N+M+2)=2/hs~2; 

A(p,N+M+3)=-2/hs~2; 

p=p+l; 

X  2nd  BC  for  Psi 

for  i=(N+M+3):(N+2*M+l) 

A(p,i-l)=l/hs~2; 
A(p,i)=-(2/hs~2-kt~2); 
A(p,i+l)=l/hs~2; 
p=p+l; 

end; 

X  3rd  BC  for  Psi 

A(p,N+M-l)=2/hs~3; 

A (p ,N+M) = (2*kl~2/hs) -4/hs~3 ; 

A(p,N+M+l)=(2/hs~3)-kt~2/hs; 

A(p,N+2*M+l)=-2/hs~2; 

A(p,N+2*M+2)=2/hs~2; 

%  Final  matrix  is  calculated 

[V,  D]=eig(A); 
D=sqrt(diag(D)); 
[Y  Il]=sort(D); 


44 


'/o  Lambda  holds  the  number  of  modes  specified  for 
'/,  use  with  later  calculations 

Lambda=Y(l rmodes) ; 

°/o  when  loop  >  1,  E_val  holds  values  for  graph 
%  of  modes  over  a  specified  frequency  range 

E_val  =  [E_val  Lambda] ; 

%  Orders  eigenvectors 

for  i  =  l:modes 

j=Ild); 

V1=[V1  V(:,j)]; 
end 

°/0  Separates  stress  potentials  for  calculations 

for  i  =  l:modes 

Phi  =  [Phi  Vl(l:N,i)]; 

Vphi  =  [Vphi  Vl(N+l:N+M+l,i)] ; 

Psi  =  [Psi  Vl(N+M+2:d,i)] ; 
end; 

'/,  Solve  for  Pressure 

Pres  =  -rf .*Phi; 

U  =  -(l/(I*w)) .*Phi; 

%  ************************************************************** 
'/,  Differencing  for  discretized  potentials 

eta2  =  Lambda. "2  -   kl~2; 
xi2  =  Lambda. "2  -  kt~2; 
zeta2  =  Lambda. ~2  -   kf~2; 

'/,   Finding  Psi_xy  (same  as  Psi_yx) 

j-i; 

for  i=l rmodes 
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Psixy(j,i)=-((Psi(j+l,i)-  (1  +  (hs-2/2).*xi2(i)).*Psi(j,i))./ 
(hs  +  (hs~3/6) .*xi2(i))); 
end 

for  j=2:M 

for  i=l :modes 

Psixy(j,i)  =  (Psi(j-l,i)-Psi(j+l,i))/(2*hs); 

end 
end 

j-M+1; 

for  i=l:modes 

Psixy(j,i)=(Psi(j-l,i)-  (1  +  (hs~2/2).*xi2(i)).*Psi(j,i))./ 

(hs  +  (hs~3/6) .*xi2(i)); 
end 

°/0  Finding  Vphi_y 

j-i; 

for  i=l:modes 

Vphiy(j,i)=-((Vphi(j+l,i)-  (1  +  (hs~2/2) .*eta2(i)) .* 

Vphi(j,i))./(hs  +  (hs~3/6).*eta2(i))); 
end 

for  j=2:M 

for  i=l rmodes 

Vphiy(j,i)  =  (Vphi(j-l,i)-Vphi(j+l,i))/(2*hs); 

end 
end 

j=M+l; 

for  i=l:modes 

Vphiy(j,i)=(Vphi(j-l,i)-  (1  +  (hs~2/2) .*eta2(i)) . * 

Vphi(j,i))./(hs  +  (hs~3/6).*eta2(i)); 
end 

X  Solving  for  Tau_yy,  Tau_xy,  and  Tau_yy 
for  i=l:modes 

Styy  =  [Styy  (-((lam  +  2*mu)*kl~2) . *Vphi( : ,i)  -  (2*mu).* 
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Psixy(: ,i))]  ; 
Stxx  =    [Stxx    ((2*mu) . *Psixy(: ,i)   -    ( (2*mu)*Lambda(i) ~2  + 

lam*kl~2) .*Vphi(: ,i))] ; 
Stxy  =    [Stxy   ( ( (2*mu*I)*Lambda(i) .*Vphiy( : , 1) )    + 

(2*Lambda(i)~2  -  kt~2)*(l/(I*Lambda(i)) ) .*Psi( : , i) )] ; 

end; 

°/,  Solving  for  u  and  v 

for   i=l:modes 

U  -    [u    (I.*Lambda(i) .*Vphi(: , i)+(l . / (I • *Lambda(i) ) ) .* 

Psixy(: ,i))]; 

v  =    [v   (Vphiy(:,i)   "  Psi(:,i))]; 
end 

"I      ************************************************************** 

%  Applying  bi-orthogonality  condition 

for  i   =   1: modes 

Fl  =    [Pres(:,i);    Styy(:,i)]; 

pr.styy  =    [pr.styy  Fl]  ; 
end 

fl   =    (Stxx'*u  -  Stxy'*v)    -  Pres'*U; 

for  j=l :modes 
for  i=l:modes 

fl(i,j)=fl(i,j)*conj(fl(iJj)); 
end 
end 

%   ************************************************************** 
°/o  Plot   of  norm  of   eigenvalues 

Yl=sqrt(Y.*conj(Y)); 

i=56:-l:l; 

figure 

plot(Yl(l:56),i,'o') 

title('Norm  of  Eigenvalues  ') 

xlabel( 'magnitude  of  lambda') 

ylabel( 'modes' ) 


47 


'/,  Determine  coefficents  for  analytical  solutions 

k=Lambda(l, l) ; 

deta2  =  k~2  -  kl~2; 

dxi2  =  k~2  -  kt~2; 

dzeta2  =  k"2  -  kf~2; 

hs=200; 

hf=10; 

h=hf +hs ; 

'/  Change  value  for  C 

C=. 00004; 

Dl=(2*k~2-kt~2)~2*tanh(sqrt(deta2)*(hs/2)); 

D2=4*k~  2*sqrt (deta2) *sqrt (dxi2) *tanh (sqrt (dxi2) * (hs/2) ) ; 

D3=(2*k~2-kt ~2) ~2*coth (sqrt (deta2) * (hs/2) ) *tanh (sqrt (dxi2) * 

(hs/2)); 

D= (D1-D2) / (D3-4*k~2*sqrt (deta2) *sqrt (dxi2) ) ; 

Bl=-(2*k~2-kt~2)/(2*I*k*sqrt(deta2)); 

B2= (cosh (sqrt (dxi2)* (hs/2)) /cosh (sqrt (deta2)* (hs/2))) *C; 

B=B1*B2; 

A2= (sinh (sqrt (dxi2) * (hs/2) ) /sinh(sqrt (det a2) * (hs/2) )  )  *D ; 

A=B1*A2; 

E1=(D1-D2)*C; 

E2=ct/(k*kt*sqrt(deta2)*(rf/rs)) ; 

E3=cosh (sqrt (dxi2) *hs/2) /sinh(sqrt (dzeta2) * (hf +hs/2) ) ; 

E=E1*E2*E3; 

y=101:-l:lJ ; 

yl=M+l:-l:l; 
y2=N:-l:l; 

Vphia  =  A.*cosh(sqrt(deta2.*yl))  +  B . *sinh(sqrt (deta2.*yl) ) ; 
Psia  =  C.*cosh(sqrt(dxi2.*yl))  +  D.*sinh(sqrt(dxi2.*yl)) ; 
Phia  =  E.*sinh(sqrt(dzeta2) .*(N-y2)) ; 

'/,  Plot  analytical  vs  discretized  solutions 

figure 

Vphia=sqrt (Vphia . *conj (Vphia) )  ; 

plot (Vphia, yl,;o' ,Vphi( : , 1) ,yl , 'x; ) 
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title( 'Analytical  vs  Discretized  for  Varphi  (solid)') 

xlabel('x') 

ylabel (' analytical    (o)      -   discretizes   (x) ' ) 

figure 

Psia  =  Psia.*I.*k; 

Psia  =  sqrt (Psia. *conj (Psia) ) ; 

plot(Psia,yl, 'o' ,Psi(: ,1) ,yl, 'x') 

title( 'Analytical  vs  Discretized  for  Psi    (solid)') 

xlabel('x') 

ylabel ( 'analytical    (o)      -   discretizes   (x)  ' ) 

figure 

Phia  =  I.*w.*Phia; 

Phia  =  sqrt (Phia. *conj (Phia)) ; 

plot (Phia, y2,'o' ,Phi( : , 1) ,y2, 'x' ) 

title( 'Analytical  vs  Discretized  for  Phi    (fluid)') 

xlabel('x') 

ylabel ( 'analytical  (o)   -  discretizes  (x)') 

yo  ************************************************************** 

'/,  end  loop 

end 
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